Sky at Night 2006 September/SAN CD 9-2006 CD-ROM 16.iso/pc/Software/Network Telescope Control/NTC-Setup.Exe/Source/ntc_server_calculus.pas

Pascal/Delphi Source File | 2006-03-24 | 36.7 KB | 1,896 lines

unit ntc_server_calculus;
{
    Copyright (C) 2004 - 2006 Andrew Sprott

    http://astronomy.crysania.co.uk
    [e-mail]

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version 2
    of the License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
}

interface

uses
    SysUtils,
    dateutils,
    math,

    ntc_server_form;

const
    null_entry=0;
    days:array[1..7] of string=(
        'monday',
        'tuesday',
        'wednesday',
        'thursday',
        'friday',
        'saturday',
        'sunday');
    light_speed=1/(3*100000);

type
    date_time=record
        hours,
        minutes,
        day,
        day_year,
        date,
        month,
        year:word;
        seconds:double;
        time_zone,
        daylight_savings:integer;
    end;

    vector=record
        x:array[1..3] of double;
    end;

    row=record
        a:array[1..3] of string[16];
    end;

    b_row=record
        b:array[1..3] of double;
    end;

    matrix=record
        q:array[1..3] of row;
        r:array[1..3] of b_row;
        binary:boolean;
        arg:double;
    end;

    tscope_calculus=class(tobject)

        { form handling }
        constructor create;

    private
        { private declarations }
        time:tdatetime;

        { utilties }
        function get_time(
            time:tdatetime)
            :date_time;

        procedure clear_time(
        period:date_time);

        { algorithms }
        function date_of_easter(
            period:date_time)
            :date_time;

    public
        { Public declarations }
        easter_sunday,
        easter_month:word;
    end;

var
    current_calculus:tscope_calculus;
    { put back }
    current_period,
    period_1980,
    period_epoc_1950,
    period_epoc_1979,
    period_epoc_1990:date_time;
    { put back }

implementation

uses
    ntc_server_observer,
    ntc_server_object;

    { -------------
        form handling
        ------------- }

constructor tscope_calculus.create;
var
    period:date_time;
begin
    inherited;
    time:=now;
    current_period:=get_time(time);
    period:=date_of_easter(current_period);
    easter_sunday:=period.date;
    easter_month:=period.month;
    with period_1980 do
        begin
            year:=1980;
            month:=1;
            date:=1;
            hours:=0;
            minutes:=0;
            seconds:=0;
        end;
    with period_epoc_1950 do
        begin
            year:=1950;
            month:=1;
            date:=1;
            hours:=0;
            minutes:=0;
            seconds:=0;
        end;
    with period_epoc_1979 do
        begin
            year:=1979;
            month:=6;
            date:=1;
            hours:=0;
            minutes:=0;
            seconds:=0;
        end;
    with period_epoc_1990 do
        begin
            year:=1990;
            month:=1;
            date:=1;
            hours:=0;
            minutes:=0;
            seconds:=0;
        end;
end;

    { --------------
        time functions
        -------------- }

procedure tscope_calculus.clear_time(
    period:date_time);
begin
    with period do
        begin
            hours:=null_entry;
            minutes:=null_entry;
            seconds:=null_entry;
            day:=null_entry;
            day_year:=null_entry;
            date:=null_entry;
            month:=null_entry;
            year:=null_entry;
            time_zone:=0;
            daylight_savings:=0;
        end;
end;

function tscope_calculus.get_time(
    time:tdatetime)
    :date_time;
var
    t:double;
begin
    result.day:=dayofweek(time);
    result.date:=dayofthemonth(time);
    result.day_year:=dayoftheyear(time);
    DecodeDate(time,result.Year,result.month,result.Day_year);
    t:=frac(time)*24;
    result.hours:=trunc(t);
    t:=frac(t)*60;
    result.seconds:=trunc(frac(t)*60);
    result.minutes:=trunc(t);
end;

    { --------------
        date of easter
        -------------- }

function tscope_calculus.date_of_easter(
    period:date_time)
    :date_time;
var
    exit_period:date_time;
    a,b,c,d,e,f,g,h,i,k,l,m,n,p:word;
    t1:word;
begin
    clear_time(exit_period);
    with period do
        begin
            a:=year mod 19;
            b:=year div 100;
            c:=year mod 100;
            d:=b div 4;
            e:=b mod 4;
            f:=(b+8) div 25;
            g:=(b-f+1) div 3;
            h:=(19*a+b-d-g+15) mod 30;
            i:=c div 4;
            k:=c mod 4;
            l:=(32+2*e+2*i-h-k) mod 7;
            m:=(a+11*h+22*i) div 451;
            t1:=(h+l-7*m+114);
            n:=t1 div 31;
            p:=t1 mod 31;
            date:=p+1;
            month:=n;
        end;
end;

    { ---------------
        day of the year
        --------------- }

function day_of_year(
    period:tdatetime)
    :word;
begin
    result:=dayoftheyear(period);
end;

    { -----------
        julian date
        ----------- }

function get_julian(
    period:date_time)
    :double;
var
    y,m,ba,bb,bc:integer;
    bd,d,t:double;
begin
    with period do
        begin
            y:=year;
            m:=month;
            t:=seconds/60;
            t:=(t+minutes)/60;
            t:=(t+hours)/24;
            d:=day+t;
            if m<3 then
                begin
                    dec(y);
                    inc(m,12);
                end;
            if (year>1582) or (month>10) or (date>=15) then
                begin
                    ba:=y div 100;
                    bb:=2-ba+ba div 4;
                end
            else
                bb:=0;
            bc:=trunc(365.25*y);
            if y<0 then
                bc:=trunc(bc-0.75);
            bd:=trunc(30.6001*(m+1));
            result:=bb+bc+bd+d+1720994.5;
        end;
end;

    { ----------------------------------
        convert julian day to calendar day
        ---------------------------------- }

function julian_to_calendar(
    julian:double)
    :date_time;
var
    bi:integer;
    bf,ba,bb,bc,bd,be,bg,d,m:double;
begin
    julian:=julian+0.5;
    bi:=trunc(julian);
    bf:=frac(julian);
    if bi>2299160 then
        begin
            ba:=trunc((bi-1867216.25)/36524.25);
            bb:=bi+1+ba-trunc(ba/4);
        end
    else
        bb:=bi;
    bc:=bb+1524;
    bd:=trunc((bc-122.1)/365.25);
    be:=trunc(365.25*bd);
    bg:=trunc((bc-be)/30.6001);
    d:=bc-be+bf-trunc(30.6001*bg);
    result.date:=trunc(d);
    d:=frac(d)*24;
    result.hours:=trunc(d);
    d:=frac(d)*60;
    result.minutes:=trunc(d);
    d:=frac(d)*60;
    result.seconds:=d;
    if bg<13.5 then
        m:=bg-1
    else
        m:=bg-13;
    if m>2.5 then
        result.year:=trunc(bd-4716)
    else
        result.year:=trunc(bd-4715);
    result.month:=trunc(m);
end;

    { -----------
        day of week
        ----------- }

function day_of_week(
    period:date_time)
    :word;
var
    jd,ba:double;
begin
    jd:=get_julian(period);
    ba:=(jd+1.5)/7;
    ba:=trunc(frac(ba)*7);
    if ba=0 then
        ba:=7;
    result:=trunc(ba);
end;

    { ---------------
        time to decimal
        --------------- }

function time_to_decimal(
    period:date_time)
    :double;
var
    t:double;
begin
    with period do
        begin
            t:=seconds/60;
            t:=(t+minutes)/60;
            t:=t+hours;
            result:=t;
        end;
end;

    { ---------------
        decimal to time
        --------------- }

function decimal_to_time(
    t:double)
    :date_time;
begin
    result.hours:=trunc(t);
    t:=frac(t)*60;
    result.minutes:=trunc(t);
    result.seconds:=trunc(frac(t)*60);
end;

    { ----------------
        local time to ut
        ---------------- }

function local_to_ut(
    period:date_time)
    :date_time;
var
    z:double;
begin
    with period do
        begin
            inc(hours,daylight_savings+time_zone);
            z:=time_to_decimal(period);
            z:=z-time_zone;
            if z>24 then
                z:=z-24;
            if z<0 then
                z:=z+24;
            result:=decimal_to_time(z);
        end;
end;

    { ----------------
        ut to local time
        ---------------- }

function ut_to_local(
    period:date_time)
    :date_time;
var
    z:double;
begin
    with period do
        begin
            z:=time_to_decimal(period);
            z:=z+time_zone;
            if z>24 then
                z:=z-24;
            if z<24 then
                z:=z+24;
            result:=decimal_to_time(z);
            result.hours:=result.hours+daylight_savings;
        end;
end;

    { -----------------------------
        ut to greenwich sidereal time
        ----------------------------- }

function ut_to_gst(
    period:date_time)
    :date_time;
var
    jd,ut,bs,bt,bt0:double;

    procedure normalise;
    begin
        while bt0>24 do
            bt0:=bt0-24;
        while bt0<24 do
            bt0:=bt0+24;
    end;

begin
    jd:=get_julian(period);
    with period do
        begin
            bs:=jd-2451545.0;
            bt:=bs/36525.0;
            bt0:=6.697374558+(2400.051336*bt)+(0.000025862*(power(bt,2)));
            normalise;
            ut:=time_to_decimal(local_to_ut(period));
            ut:=ut*1.002737090;
            bt0:=bt0+ut;
            normalise;
            result:=decimal_to_time(bt0);
        end;
end;

    { -----------------------------
        greenwich sidereal time to ut
        ----------------------------- }

function gst_to_ut(
    period:date_time)
    :date_time;
var
    jd,ut,bs,bt,bt0:double;

    procedure normalise;
    begin
        while bt0>24 do
            bt0:=bt0-24;
        while bt0<24 do
            bt0:=bt0+24;
    end;

begin
    jd:=get_julian(period);
    with period do
        begin
            bs:=jd-2451545.0;
            bt:=bs/36525.0;
            bt0:=6.697374558+(2400.051336*bt)+(0.000025862*(power(bt,2)));
            normalise;
            ut:=time_to_decimal(period);
            normalise;
            ut:=ut*0.9972695663;
            result:=decimal_to_time(ut);
        end;
end;

    { ----------
        gst to lst
        ---------- }

function gst_to_lst(
    period:date_time)
    :date_time;
var
    gst,l:double;
begin
    gst:=time_to_decimal(period);
    with scope_observer do
        begin
            l:=longitude/15;
            if longitude>0 then
                gst:=gst-l
            else if longitude<0 then
                gst:=gst+l;
            if gst>24 then
                gst:=gst-24;
            if gst<0 then
                gst:=gst+24;
        end;
    result:=decimal_to_time(gst);
end;

    { ----------
        lst to gst
        ---------- }

function lst_to_gst_decimal(
    lst:double)
    :double;
var
    l:double;
begin
    with scope_observer do
        begin
            l:=longitude/15;
            if longitude>0 then
                lst:=lst+l
            else if longitude<0 then
                lst:=lst-l;
            if lst>24 then
                lst:=lst-24;
            if lst<0 then
                lst:=lst+24;
        end;
    result:=lst;
end;

function lst_to_gst(
    period:date_time)
    :date_time;
var
    lst,gst:double;
begin
    lst:=time_to_decimal(period);
    gst:=lst_to_gst_decimal(lst);
    result:=decimal_to_time(gst);
end;

    { ------------------
        degrees to decimal
        ------------------ }

function degrees_to_decimal(
    coord:dms)
    :double;
var
    d:double;
begin
    with coord do
        begin
            d:=seconds/60;
            d:=(d+minutes)/60;
            d:=d+degrees;
            result:=d;
        end;
end;

    { ------------------
        decimal to degrees
        ------------------ }

function decimal_to_degrees(
    d:double)
    :dms;
var
    t:double;
begin
    result.degrees:=trunc(d);
    t:=frac(d)*60;
    result.minutes:=trunc(t);
    result.seconds:=trunc(frac(t)*60);
    result.minutes:=trunc(t);
end;

    { -------------
        ra to decimal
        ------------- }

function ra_to_decimal(
    coord:hms)
    :double;
var
    t:double;
begin
    with coord do
        begin
            t:=seconds/60;
            t:=(t+minutes)/60;
            t:=t+hours;
            result:=t;
        end;
end;

    { -------------
        decimal to ra
        ------------- }

function decimal_to_ra(
    t:double)
    :hms;
var
    d:double;
begin
    result.hours:=trunc(t);
    d:=frac(t)*60;
    result.minutes:=trunc(d);
    result.seconds:=trunc(frac(d)*60);
    result.minutes:=trunc(d);
end;

    { ----------------------
        right ascension to lst
        ---------------------- }

function ra_to_lst(
    period:date_time;
    ra:hms)
    :date_time;
var
    gst,lst:date_time;
    ra_d,t,lst_d:double;
begin
    gst:=ut_to_gst(period);
    lst:=gst_to_lst(gst);
    lst_d:=time_to_decimal(lst);
    ra_d:=ra_to_decimal(ra);
    t:=lst_d-ra_d;
    if t<0 then
        t:=t+24;
    result:=decimal_to_time(t);
end;

    { ----------------------
        lst to right ascension
        ---------------------- }

function lst_to_ra(
    ra:hms;
    period:date_time)
    :hms;
var
    lst:date_time;
    t,lst_d:double;
begin
    lst:=gst_to_lst(ut_to_gst(period));
    with period do
        begin
            hours:=trunc(ra_to_decimal(ra));
            lst_d:=time_to_decimal(lst);
            t:=lst_d-hours;
            if t<0 then
                t:=t+24;
        end;
    result:=decimal_to_ra(t);
end;

    { ------------------
        ecliptic obliquity
        ------------------ }

function ecliptic_obliquity
    :double;
var
    jd,bt,ao:double;
begin
    jd:=get_julian(period_1980);
    jd:=jd-2451545.0;
    bt:=jd/36525.0;
    ao:=46.815*bt+0.0006*power(bt,2)-0.00181*power(bt,3);
    ao:=ao/3600;
    result:=23.439292-ao;
end;

    { ------------------------------------
        match degree with signs of quadrants
        ------------------------------------ }

procedure match_xy_with_degree(
            x,
            y:double;
    var deg:double);
begin
    if x>=0 then
        begin
            if y>=0 then
                begin
                    if deg<-360 then
                        deg:=deg+360
                    else if deg<-180 then
                        deg:=deg+180
                    else if deg>360 then
                        deg:=deg-360
                    else if deg>180 then
                        deg:=deg-180;
                end
            else
                begin
                    if deg<-270 then
                        deg:=deg+360
                    else if deg<-90 then
                        deg:=deg+180
                    else if deg>270 then
                        deg:=deg-360
                    else if deg>90 then
                        deg:=deg-180;
                end;
        end
    else
        begin
            if y<0 then
                begin
                    if deg<-180 then
                        deg:=deg+360
                    else if deg<-360 then
                        deg:=deg+180
                    else if deg>180 then
                        deg:=deg-360
                    else if deg>360 then
                        deg:=deg-180;
                end
            else
                begin
                    if deg<-90 then
                        deg:=deg+360
                    else if deg<-270 then
                        deg:=deg+180
                    else if deg>90 then
                        deg:=deg-360
                    else if deg>270 then
                        deg:=deg-180;
                end;
        end;
end;

    { --------
        matrices
        -------- }

function setup_v(
    u,
    v:double)
    :vector;
begin
    result.x[1]:=cos(u)*cos(v);
    result.x[2]:=sin(u)*cos(v);
    result.x[3]:=sin(v);
end;

function get_ra_dec_vector(
    ra:hms;
    dec:dms)
    :vector;
var
    u,v:double;
begin
    u:=ra_to_decimal(ra);
    v:=degrees_to_decimal(dec);
    result:=setup_v(u,v);
end;

function get_degrees_vector(
    u_d:dms;
    v_d:dms)
    :vector;
var
    u,v:double;
begin
    u:=degrees_to_decimal(u_d);
    v:=degrees_to_decimal(v_d);
    result:=setup_v(u,v);
end;

procedure process_matrix(
            w:vector;
            m:matrix;
    var v:vector;
            arg:double);
var
    i:integer;

    function parse(
        s:string)
        :double;

        function check_function(
            i:integer)
            :double;
        var
            t:string;
        begin
            result:=0;
            t:=copy(s,i+5,3);
            if i=1 then
                begin
                    if t='cos' then
                        result:=cos(arg)
                    else if t='sin' then
                        result:=sin(arg);
                end
            else
                begin
                    if t='cos' then
                        result:=-cos(arg)
                    else if t='sin' then
                        result:=-sin(arg);
                end;
        end;

    begin
        if s[1]='-' then
            begin
                if (s[2]>='0') and (s[2]<='9') then
                    result:=strtofloatdef(s,0)
                else
                    result:=check_function(2);
            end
        else if (s[1]>='0') and (s[1]<='9') then
            result:=strtofloatdef(s,0)
        else
            result:=check_function(1);
    end;

begin
    with m do
         if not binary then
        begin
            for i:=1 to 3 do
                 with q[i],v do
                w.x[i]:=parse(a[1])*x[1]+parse(a[2])*x[2]+parse(a[3])*x[3];
        end
    else
        begin
            for i:=1 to 3 do
                 with r[i],v do
                w.x[i]:=b[1]*x[1]+b[2]*x[2]+b[3]*x[3];
        end;
end;

procedure get_vector_ra_dec(
    w:vector;
    ra:hms;
    dec:dms);
var
    u,v:double;
begin
    with w do
        begin
            u:=arctan(x[2]/x[1]);
            v:=arcsin(x[3]);
        end;
    ra:=decimal_to_ra(u);
    dec:=decimal_to_degrees(v);
end;

procedure get_vector_degrees(
    w:vector;
    lat,
    long:dms);
var
    u,v:double;
begin
    with w do
        begin
            u:=arctan(x[2]/x[1]);
            v:=arcsin(x[3]);
        end;
    lat:=decimal_to_degrees(u);
    long:=decimal_to_degrees(v);
end;

    { --------------------------
        local position and horizon
        -------------------------- }

procedure horizon_to_localised(
            alt,
            az:dms;
    var hours:hms;
    var    dec:dms);
var
    v,w:vector;
    m:matrix;
    latitude:double;
begin
    v:=get_degrees_vector(alt,az);
    latitude:=scope_observer.latitude;
    with m do
        begin
            arg:=latitude;
            binary:=false;
            q[1].a[1]:='-sin';
            q[1].a[2]:='';
            q[1].a[3]:='cos';
            q[2].a[1]:='';
            q[2].a[2]:='-1';
            q[2].a[3]:='';
            q[3].a[1]:='cos';
            q[3].a[2]:='';
            q[3].a[3]:='sin';
            process_matrix(w,m,v,arg);
        end;
    get_vector_ra_dec(w,hours,dec);
end;

procedure localised_to_horizon(
            hours:hms;
            dec:dms;
    var alt,
            az:dms);
var
    v,w:vector;
    m:matrix;
    latitude:double;
begin
    v:=get_ra_dec_vector(hours,dec);
    latitude:=scope_observer.latitude;
    with m do
        begin
            arg:=latitude;
            binary:=false;
            q[1].a[1]:='-sin';
            q[1].a[2]:='';
            q[1].a[3]:='cos';
            q[2].a[1]:='';
            q[2].a[2]:='-1';
            q[2].a[3]:='';
            q[3].a[1]:='cos';
            q[3].a[2]:='';
            q[3].a[3]:='sin';
            process_matrix(w,m,v,arg);
        end;
    get_vector_degrees(w,alt,az);
end;

    { ------------------------
        localised and equatorial
        ------------------------ }

procedure equatorial_to_localised(
            ra:hms;
            dec:dms;
    var l_ra:hms;
    var l_dec:dms);
var
    v,w:vector;
    m:matrix;
    lst:double;
begin
    v:=get_ra_dec_vector(ra,dec);
    lst:=time_to_decimal(ra_to_lst(current_period,ra));
    with m do
        begin
            arg:=lst;
            binary:=false;
            q[1].a[1]:='cos';
            q[1].a[2]:='sin';
            q[1].a[3]:='';
            q[2].a[1]:='sin';
            q[2].a[2]:='-cos';
            q[2].a[3]:='';
            q[3].a[1]:='';
            q[3].a[2]:='';
            q[3].a[3]:='1';
            process_matrix(w,m,v,arg);
        end;
    get_vector_ra_dec(w,l_ra,l_dec);
end;

procedure localised_to_equatorial(
            l_ra:hms;
            l_dec:dms;
    var    ra:hms;
    var dec:dms);
var
    v,w:vector;
    m:matrix;
    lst:double;
begin
    v:=get_ra_dec_vector(l_ra,l_dec);
    lst:=time_to_decimal(ra_to_lst(current_period,ra));
    with m do
        begin
            arg:=lst;
            binary:=false;
            q[1].a[1]:='cos';
            q[1].a[2]:='sin';
            q[1].a[3]:='';
            q[2].a[1]:='sin';
            q[2].a[2]:='-cos';
            q[2].a[3]:='';
            q[3].a[1]:='';
            q[3].a[2]:='';
            q[3].a[3]:='1';
            process_matrix(w,m,v,arg);
        end;
    get_vector_ra_dec(w,ra,dec);
end;

    { -----------------------
        ecliptic and equatorial
        ----------------------- }

procedure equatorial_to_ecliptic(
            ra:hms;
            dec:dms;
    var e_long,
            e_lat:dms);
var
    v,w:vector;
    m:matrix;
    obliquity:double;
begin
    v:=get_ra_dec_vector(ra,dec);
    obliquity:=ecliptic_obliquity;
    with m do
        begin
            arg:=obliquity;
            binary:=false;
            q[1].a[1]:='1';
            q[1].a[2]:='';
            q[1].a[3]:='';
            q[2].a[1]:='';
            q[2].a[2]:='cos';
            q[2].a[3]:='sin';
            q[3].a[1]:='';
            q[3].a[2]:='-sin';
            q[3].a[3]:='cos';
            process_matrix(w,m,v,arg);
        end;
    get_vector_degrees(w,e_long,e_lat);
end;

procedure ecliptic_to_equatorial(
            e_long,
            e_lat:dms;
    var ra:hms;
    var dec:dms);
var
    v,w:vector;
    m:matrix;
    obliquity:double;
begin
    v:=get_degrees_vector(e_long,e_lat);
    obliquity:=ecliptic_obliquity;
    with m do
        begin
            arg:=obliquity;
            binary:=false;
            q[1].a[1]:='1';
            q[1].a[2]:='';
            q[1].a[3]:='';
            q[2].a[1]:='';
            q[2].a[2]:='cos';
            q[2].a[3]:='-sin';
            q[3].a[1]:='';
            q[3].a[2]:='sin';
            q[3].a[3]:='cos';
            process_matrix(w,m,v,arg);
        end;
    get_vector_ra_dec(w,ra,dec);
end;

    { ----------------------
        glactic and equatorial
        ---------------------- }

procedure equatorial_to_galactic(
            ra:hms;
            dec:dms;
    var g_long,
            g_lat:dms);
var
    v,w:vector;
    m:matrix;
begin
    v:=get_ra_dec_vector(ra,dec);
    with m do
        begin
            arg:=0;
            binary:=false;
            q[1].a[1]:='-0.0669887';
            q[1].a[2]:='-0.8727558';
            q[1].a[3]:='-0.4835389';
            q[2].a[1]:='0.4927285';
            q[2].a[2]:='-0.4503470';
            q[2].a[3]:='0.7445846';
            q[3].a[1]:='-0.8676008';
            q[3].a[2]:='-0.1883746';
            q[3].a[3]:='0.4601998';
            process_matrix(w,m,v,arg);
        end;
    get_vector_degrees(w,g_long,g_lat);
end;

procedure galactic_to_equatorial(
            g_long,
            g_lat:dms;
    var ra:hms;
    var dec:dms);
var
    v,w:vector;
    m:matrix;
begin
    v:=get_degrees_vector(g_long,g_lat);
    with m do
        begin
            arg:=0;
            binary:=false;
            q[1].a[1]:='-0.0669887';
            q[1].a[2]:='-0.8727558';
            q[1].a[3]:='-0.8676008';
            q[2].a[1]:='-0.8727558';
            q[2].a[2]:='-0.4503470';
            q[2].a[3]:='-0.1883746';
            q[3].a[1]:='-0.4835389';
            q[3].a[2]:='0.7445846';
            q[3].a[3]:='0.4601998';
            process_matrix(w,m,v,arg);
        end;
    get_vector_ra_dec(w,ra,dec);
end;

    { ---------------------------------------------
        find angle between two equatorial coordinates
        --------------------------------------------- }

function find_angle(
    ra_1,
    ra_2:hms;
    dec_1,
    dec_2:dms)
    :dms;
var
    ra_1_d,ra_2_d,
    dec_1_d,dec_2_d,
    t,d:double;
begin
    ra_1_d:=ra_to_decimal(ra_1);
    dec_1_d:=degrees_to_decimal(dec_1);
    ra_2_d:=ra_to_decimal(ra_2);
    dec_2_d:=degrees_to_decimal(dec_2);
    t:=(ra_1_d-ra_2_d)*15;
    d:=sin(dec_1_d)*sin(dec_2_d)+cos(dec_1_d)*cos(dec_2_d)*cos(t);
    d:=arccos(d);
    result:=decimal_to_degrees(d);
end;

    { --------------------------------------------
        find rising and setting times for coordinate
        -------------------------------------------- }

function rise_and_set(
            ra:hms;
            dec:dms;
    var rising_d,
            setting_d:double;
    var rising_az,
            setting_az:dms;
    var ang:double)
            :boolean;
var
    ra_d,
    dec_d,
    a1,a2,
    h,u,x,y:double;
begin
    ra_d:=ra_to_decimal(ra);
    dec_d:=degrees_to_decimal(dec);
    a1:=(sin(dec_d)/cos(scope_observer.latitude));
    if (a1<1) and (a1>-1) then
        begin
            a1:=arccos(a1);
            result:=true;
            a2:=360-a1;
            u:=arccos(sin(scope_observer.latitude)/cos(a2));
            x:=34/60;
            { calculate refraction }
            ang:=arcsin(tan(x)/tan(u));
            a1:=a1-ang;
            a2:=a2+ang;
            rising_az:=decimal_to_degrees(a1);
            setting_az:=decimal_to_degrees(a2);
            h:=(1/15)*arccos(-tan(scope_observer.latitude)*tan(dec_d));
            a1:=24+ra_d+h;
            if a1>24 then
                a1:=a1-24;
            a2:=ra_d+h;
            if a2>24 then
                a2:=a2-24;
            { calculate refraction }
            y:=arcsin(sin(x)/sin(u));
            ang:=(240/y)/cos(dec_d);
            ang:=ang/3600;
            a1:=a1-ang;
            a2:=a2+ang;
            rising_d:=lst_to_gst_decimal(a1);
            setting_d:=lst_to_gst_decimal(a2);
        end
    else
        result:=false;
end;

    { ----------
        precession
        ---------- }

procedure precession(
            ra_1:hms;
            dec_1:dms;
    var ra_2:hms;
    var dec_2:dms);
var
    jd1,p1,p2,p3,
    ra_1_d,dec_1_d,
    u2,v2:double;
var
    v,w,s:vector;
    m:matrix;

    procedure precession_variables;
    begin
        p1:=0.6406161*jd1+0.0000839*power(jd1,2)+0.0000050*power(jd1,3);
        p2:=0.6406161*jd1+0.0003041*power(jd1,2)+0.0000051*power(jd1,3);
        p3:=0.5567530*jd1-0.0001185*power(jd1,2)-0.0000116*power(jd1,3);
    end;

begin
    jd1:=get_julian(period_epoc_1950);
    jd1:=(jd1-2451545)/36525;
    precession_variables;
    with m do
        begin
            binary:=true;
            arg:=0;
            { cx*ct*cz-sx*sz }
            r[1].b[1]:=cos(p1)*cos(p3)*cos(p2)-sin(p1)*sin(p2);
            { cx*ct*sz+sx*cz }
            r[1].b[2]:=cos(p1)*cos(p3)*sin(p2)+sin(p1)*cos(p2);
            { cx*st }
            r[1].b[3]:=cos(p1)*sin(p3);
            { -sx*ct*cz-cx*sz }
            r[2].b[1]:=-sin(p1)*cos(p3)*cos(p2)-cos(p1)*sin(p2);
            { -sx*ct*sz+cx*cz }
            r[2].b[2]:=-sin(p1)*cos(p3)*sin(p2)+cos(p1)*cos(p2);
            { -sx*st }
            r[2].b[3]:=-sin(p1)*sin(p3);
            { -st*cz }
            r[3].b[1]:=-sin(p3)*cos(p2);
            { -st*sz }
            r[3].b[2]:=-sin(p3)*sin(p2);
            { ct }
            r[3].b[3]:=cos(p3);
        end;
    ra_1_d:=ra_to_decimal(ra_1);
    dec_1_d:=degrees_to_decimal(dec_1);
    v:=setup_v(ra_1_d,dec_1_d);
    process_matrix(s,m,v,0);
    jd1:=get_julian(period_epoc_1979);
    jd1:=(jd1-2451545)/36525;
    precession_variables;
    with m do
        begin
            { cx*ct*cz-sx*sz }
            r[1].b[1]:=cos(p1)*cos(p3)*cos(p2)-sin(p1)*sin(p2);
            { -sx*ct*cz-cx*sz }
            r[1].b[2]:=-sin(p1)*cos(p3)*cos(p2)-cos(p1)*sin(p2);
            { -st*cz }
            r[1].b[3]:=-sin(p3)*cos(p2);
            { cx*ct*sz+sx*cz }
            r[2].b[1]:=cos(p1)*cos(p3)*sin(p2)+sin(p1)*cos(p2);
            { -sx*ct*sz+cx*cz }
            r[2].b[2]:=-sin(p1)*cos(p3)*sin(p2)+cos(p1)*cos(p2);
            { -st*sz }
            r[2].b[3]:=-sin(p3)*sin(p2);
            { cx*st }
            r[3].b[1]:=cos(p1)*sin(p3);
            { -sx*st }
            r[3].b[2]:=-sin(p1)*sin(p3);
            { ct }
            r[3].b[3]:=cos(p3);
        end;
    process_matrix(w,m,s,0);
    with w do
        begin
            u2:=arctan(x[2]/x[1]);
            v2:=arcsin(x[3]);
            match_xy_with_degree(x[1],x[2],u2);
        end;
    ra_2:=decimal_to_ra(u2);
    dec_2:=decimal_to_degrees(v2);
end;

    { --------
        nutation
        -------- }

procedure nutation(
    period:date_time;
    var e_long_a,
            e_ob_a:double);
var
    jd,a,b,l,ma:double;

    procedure keep_in_range(
        var f:double);
    begin
        while f>360 do
            f:=f-360;
        while f<0 do
            f:=f+360;
    end;

begin
    jd:=get_julian(period);
    jd:=(jd-2415020)/36525;
    a:=100.002136*jd;
    l:=279.6967+360*(a-trunc(a));
    keep_in_range(l);
    b:=5.372617*jd;
    ma:=259.1833-360*(b-trunc(b));
    keep_in_range(ma);
    e_long_a:=-17.2*sin(ma)-1.3*sin(2*l);
    e_ob_a:=9.2*cos(ma)+0.5*cos(2*l);
end;

    { ----------
        aberration
        ---------- }

procedure aberration(
    var e_long,
            e_lat:dms;
            s_e_long:dms);
var
    e_long_d,e_lat_d,s_e_long_d:double;
begin
    e_long_d:=degrees_to_decimal(e_long);
    e_lat_d:=degrees_to_decimal(e_lat);
    s_e_long_d:=degrees_to_decimal(s_e_long);
    e_long_d:=e_long_d+(-20.5*cos(s_e_long_d-e_long_d)/cos(e_lat_d))/3600;
    e_lat_d:=e_lat_d+(-20.5*sin(s_e_long_d-e_long_d)*sin(e_lat_d))/3600;
    e_long:=decimal_to_degrees(e_long_d);
    e_lat:=decimal_to_degrees(e_lat_d);
end;

    { ----------
        refraction
        ---------- }

procedure refraction(
            alt:dms;
            p,
            t:double;
    var a_alt:dms);
var
    z,alt_d:double;
begin
    alt_d:=degrees_to_decimal(alt);
    if alt_d>15 then
        begin
            z:=90-alt_d;
            alt_d:=alt_d+(0.00452*p*tan(z)/(273+t));
        end
    else
         alt_d:=alt_d+
        (p*(0.1594+0.0196*alt_d+0.00002*power(alt_d,2)))/
        ((273+t)*(1+0.505*alt_d+0.0845*power(alt_d,2)));
    a_alt:=decimal_to_degrees(alt_d);
end;

    { -------------------
        geocentric parallax
        ------------------- }

procedure geocentric_parallax(
            lat:dms;
            h:double;
    var p_sin_l,
            p_cos_l:double);
var
    u,lat_d:double;
begin
    lat_d:=degrees_to_decimal(lat);
    u:=arctan(0.996647*tan(lat_d));
    h:=(h/6378140);
    p_sin_l:=0.996647*sin(u);
    p_cos_l:=cos(u)+h*cos(lat_d);
end;

    { ------------------------------------------
        correct for geocentric parallax for object
        ------------------------------------------ }

procedure correct_for_parallax(
            period:date_time;
            au,
            height:double;
            lat:dms;
            ra:hms;
            dec:dms;
    var a_ra:hms;
    var a_dec:dms);
var
    lst:date_time;
    pie,h,a1,a2,
    ra_d,ra_dd,dec_d,dec_dd,p_cos_l,p_sin_l:double;
begin
    pie:=(8.794/au)/3600/15;
    lst:=gst_to_lst(ut_to_gst(period));
    ra_d:=ra_to_decimal(ra);
    dec_d:=degrees_to_decimal(dec);
    geocentric_parallax(lat,height,p_sin_l,p_cos_l);
    h:=ra_to_decimal(lst_to_ra(ra,lst));
    a1:=(pie*sin(h)*p_cos_l)/cos(dec_d);
    ra_dd:=ra_d-a1;
    a2:=(pie*(p_sin_l*cos(dec_d)-p_cos_l*cos(h)*sin(dec_d)));
    dec_dd:=dec_d-a2;
    a_ra:=decimal_to_ra(ra_dd);
    a_dec:=decimal_to_degrees(dec_dd);
end;

    { ===========================================
        planets, asteroids, binary stars and comets
        =========================================== }



    { ==============================
        heliographic coordinates [sun]
        ============================== }

    { -----------------------
        calculate suns position
        ----------------------- }

procedure suns_position(
            period:date_time;
    var ra:hms;
    var dec:dms;
    var distance,
            size,
            e_long_d:double);
var
    jd,jd_1990,n,m,e_long_e1990_d,e_long_p_d,ecc_d,ecc_r,
    ec,a,e,ae,v,f,sapc:double;
    done:boolean;
begin
    jd_1990:=get_julian(period_epoc_1990);
    jd:=get_julian(period)-jd_1990;
    n:=(360/365.242191)*trunc(jd);
    while n<0 do
        n:=n+360;
    while n>360 do
        n:=n-360;
    e_long_e1990_d:=279.6966778+36000.76892+0.0003025*power(jd_1990,2);
    e_long_p_d:=281.2208444+1.719175*jd_1990+0.000452778*power(jd_1990,2);
    ecc_d:=0.01675104-0.0000418*jd_1990-0.000000126*power(jd_1990,2);
    m:=n+e_long_e1990_d-e_long_p_d;
    if ecc_d<=0.1 then
        begin
            ecc_r:=ecc_d*(pi/180);
            e:=m;
            done:=false;
            while not done do
                begin
                    a:=e-ecc_r*sin(e)-m;
                    done:=a<=0.000001;
                    if not done then
                        begin
                            ae:=a/(1-ecc_r*cos(e));
                            e:=ae;
                        end;
                end;
            v:=sqrt((1+ecc_r)/(1-ecc_r))*tan(e/2);
            v:=(arctan(v)*2)*(180/pi);
            e_long_d:=v+e_long_e1990_d;
            if e_long_d<0 then
                e_long_d:=e_long_d+360
            else if e_long_d>360 then
                e_long_d:=e_long_d-360;
        end
    else
        begin
            if m<0 then
                m:=m+360;
            ec:=360/pi*ecc_d*sin(m);
            v:=m+ec;
            e_long_d:=n+ec+e_long_e1990_d;
            if e_long_d>360 then
                e_long_d:=e_long_d-360;
        end;
    f:=(1+ecc_d*cos(v))/(1-power(ecc_d,2));
    distance:=(1.495985/f)*power(10,8);
    size:=f*0.533128;
    sapc:=distance/light_speed;
    sapc:=360/365.242191*(sapc/60/60/24);
    e_long_d:=e_long_d-sapc;
    ecliptic_to_equatorial(
        decimal_to_degrees(e_long_d),
        decimal_to_degrees(0),
        ra,dec);
end;

    { -----------------------
        find coordinates of sun
        ----------------------- }

procedure find_sunspot(
            period:date_time;
            pos_angle,
            displacement:dms;
    var h_long_c,
            h_lat_c:double;
    var h_long,
            h_lat:dms);
var
    jd,t,an_long,x,y,a,m,ad1,ad2,o,
    p_axis,pos_ang_d,dis_d,rad_d,p,size,
    h_lat_d,h_long_d,distance,gd:double;
    ra:hms;
    dec,tp:dms;
begin
    jd:=get_julian(period);
    t:=(jd-2415020)/36525;
    a:=84*t/60;
    tp.degrees:=74;
    tp.minutes:=22;
    tp.seconds:=0;
    suns_position(period,ra,dec,distance,size,gd);
    an_long:=degrees_to_decimal(tp)+a;
    y:=sin(an_long-gd)*cos(7.25);
    x:=-cos(an_long-gd);
    a:=arctan(y/x);
    match_xy_with_degree(x,y,a);
    m:=(360/25.38)*(jd-2398220);
    while m>360 do
        m:=m-360;
    m:=360-m;
    a:=a+m;
    if a>360 then
        a:=a-360;
    h_long_d:=a;
    h_lat_d:=arcsin(sin(gd-an_long)*sin(7.25));
    h_long_c:=h_long_d;
    h_lat_c:=h_lat_d;
    o:=ecliptic_obliquity;
    ad1:=arctan(-cos(gd)*tan(o));
    ad2:=arctan(-cos(an_long-gd)*tan(7.25));
    p_axis:=ad1+ad2;
    dis_d:=degrees_to_decimal(displacement);
    rad_d:=size/2;
    p:=arcsin(dis_d/rad_d)-dis_d;
    pos_ang_d:=degrees_to_decimal(pos_angle);
    h_lat_d:=
        arcsin((sin(h_lat_d)*cos(p)+cos(h_lat_d)*sin(p)*cos(p_axis-pos_ang_d)));
    a:=arcsin((sin(p)*sin(p_axis-pos_ang_d)/(cos(h_lat_d))));
    h_long_d:=h_long_d+a;
    if h_long_d>360 then
        h_long_d:=h_long_d-360;
    h_lat:=decimal_to_degrees(h_lat_d);
    h_long:=decimal_to_degrees(h_long_d);
end;

    { ------------------
        sunrise and sunset
        ------------------ }

procedure sunrise_and_sunset(
            period:date_time;
            longitude,
            latitude:dms;
            height:double;
    var rises,
            sets:date_time);
var
    midnight:date_time;
    ra,ra_b,ra_a,a_ra:hms;
    dec,dec_b,a_dec,dec_a,blank:dms;
    distance,size,ang,ang_1,ang_2,a_lat,c,x,a_dec_d,
    gst_r_1,gst_r_2,gst_s_1,gst_s_2,gst_r,gst_s,
    gd,gd_b,gd_a,long_d,gst_0,gst_0_ut:double;
    zero_lat:dms;
begin
    midnight:=period;
    with midnight do
        begin
            hours:=0;
            minutes:=0;
            seconds:=0;
        end;
    suns_position(period,ra,dec,distance,size,gd);
    gd_b:=gd-0.985647;
    with zero_lat do
        begin
            degrees:=0;
            minutes:=0;
            seconds:=0;
        end;
    ecliptic_to_equatorial(decimal_to_degrees(gd_b),zero_lat,ra_b,dec_b);
    rise_and_set(ra_b,dec_b,gst_r_1,gst_s_1,blank,blank,ang_1);
    gd_a:=gd+0.985647;
    ecliptic_to_equatorial(decimal_to_degrees(gd_a),zero_lat,ra_a,dec_a);
    rise_and_set(ra_a,dec_a,gst_r_2,gst_s_2,blank,blank,ang_2);
    if gd_b>gd_a then
        begin
            gst_r_2:=gst_r_2+24;
            gst_s_2:=gst_s_2+24;
        end;
    gst_0:=time_to_decimal(ut_to_gst(period));
    long_d:=degrees_to_decimal(longitude);
    long_d:=long_d/15*1.002738;
    gst_0_ut:=gst_0-long_d;
    if gst_0_ut<0 then
        gst_0_ut:=gst_0_ut+24;
    if gd_b<gst_0_ut then
        begin
            gst_r_1:=gst_r_1+24;
            gst_s_1:=gst_s_1+24;
            gst_r_2:=gst_r_2+24;
            gst_s_2:=gst_s_2+24;
        end;
    gst_r:=(24.07*gst_r_1-gst_0*(gst_r_2-gst_r_1))/(24.07+gst_r_1-gst_r_2);
    gst_s:=(24.07*gst_s_1-gst_0*(gst_s_2-gst_s_1))/(24.07+gst_s_1-gst_s_2);
    correct_for_parallax(period,distance,height,latitude,ra,dec,a_ra,a_dec);
    a_lat:=(degrees_to_decimal(dec_b)+degrees_to_decimal(dec_a))/2;
    ang:=(ang_1+ang_2)/2;
    a_dec_d:=degrees_to_decimal(a_dec);
    x:=(size/2)/(-a_dec_d+ang);
    c:=-(a_lat+ang+x)/3600;
    gst_r:=gst_r-c;
    gst_s:=gst_s+c;
    rises:=midnight;
    with rises do
        begin
            rises:=decimal_to_time(gst_r);
            rises:=gst_to_ut(rises);
            rises:=ut_to_local(rises);
        end;
    sets:=midnight;
    with sets do
        begin
            sets:=decimal_to_time(gst_s);
            sets:=gst_to_ut(sets);
            sets:=ut_to_local(sets);
        end;
end;

procedure twilight(
            period:date_time;
    var evening,
            morning:date_time;
    var all_night:boolean);
var
    ra:hms;
    dec:dms;
    rises,sets:date_time;
    h1,h2,t,r,s,
    dec_d,distance,size,e_long_d:double;
begin
    with scope_observer do
        begin
            with period do
                begin
                    hours:=12;
                    minutes:=0;
                    seconds:=0;
                end;
            suns_position(period,ra,dec,distance,size,e_long_d);
            dec_d:=degrees_to_decimal(dec);
            h1:=arccos(-tan(latitude)*tan(dec_d));
            h2:=(cos(108)-sin(latitude)*sin(dec_d))/(cos(latitude)*cos(dec_d));
            if (h2<1) and (h2>-1) then
                begin
                    all_night:=false;
                    h2:=arccos(h2);
                t:=((h2-h1)/15)*0.9973;
                    sunrise_and_sunset(period,
                        decimal_to_degrees(longitude),
                        decimal_to_degrees(latitude),
                        height,rises,sets);
                    r:=time_to_decimal(rises);
                    r:=r-t;
                    morning:=decimal_to_time(r);
                    s:=time_to_decimal(sets);
                    s:=s+t;
                    evening:=decimal_to_time(s);
                end
            else
                begin
                    morning:=period;
                    evening:=period;
                    all_night:=true;
                end;
        end;
end;

    { ------------------------------------
        calculate carrington rotation number
        ------------------------------------ }

function get_crn(
    period:date_time)
    :integer;
var
    jd:double;
begin
    jd:=get_julian(period);
    result:=round(1690+((jd-2444235.34)/272753));
end;

    { ----------------
        equation of time
        ---------------- }

function equation_of_time(
    period:date_time)
    :hms;
var
    ra:hms;
    dec:dms;
    midday:date_time;
    distance,size,e_long_d,ut:double;
begin
    suns_position(period,ra,dec,distance,size,e_long_d);
    with period do
        begin
            hours:=ra.hours;
            minutes:=ra.minutes;
            seconds:=ra.seconds;
        end;
    midday:=period;
    with midday do
        begin
            hours:=12;
            minutes:=0;
            seconds:=0;
        end;
    ut:=time_to_decimal(gst_to_ut(midday))-time_to_decimal(gst_to_ut(period));
    midday:=decimal_to_time(ut);
    with midday do
        begin
            result.hours:=hours;
            result.minutes:=minutes;
            result.seconds:=trunc(seconds);
        end;
end;

    { ----------------
        solar elongation
        ---------------- }

function solar_elongation(
    period:date_time;
    ra_p:hms;
    dec_p:dms)
    :double;
var
    ra:hms;
    dec:dms;
    distance,size,e_long_d,ra_d,dec_d,ra_p_d,dec_p_d:double;
begin
    suns_position(period,ra,dec,distance,size,e_long_d);
    ra_d:=ra_to_decimal(ra);
    dec_d:=degrees_to_decimal(dec);
    ra_p_d:=ra_to_decimal(ra_p);
    dec_p_d:=degrees_to_decimal(dec_p);
    ra_p_d:=ra_p_d*15;
    result:=arccos(
        (sin(dec_p_d)*sin(dec_d)+cos(ra_p_d-ra_d)*cos(dec_p_d)*cos(dec_d)));
end;

    { ================================
        selenographic coordinates [moon]
        ================================ }

    { --------------------------------------------
        correct for geocentric parallax for the moon
        -------------------------------------------- }

procedure correct_for_parallax_moon(
            period:date_time;
            height:double;
            moon_p,
            lat:dms;
            ra:hms;
            dec:dms;
    var a_ra:hms;
    var a_dec:dms);
var
    lst:date_time;
    a,r,h,hh,ra_d,ra_dd,dec_d,dec_dd,moon_p_d,
    p_sin_l,p_cos_l:double;
begin
    ra_d:=ra_to_decimal(ra);
    dec_d:=degrees_to_decimal(dec);
    lst:=gst_to_lst(ut_to_gst(period));
    h:=ra_to_decimal(lst_to_ra(ra,lst));
    geocentric_parallax(lat,height,p_sin_l,p_cos_l);
    moon_p_d:=degrees_to_decimal(moon_p);
    r:=(1/sin(moon_p_d));
    a:=arctan((p_cos_l*sin(h))/(r*cos(dec_d)-p_cos_l*cos(h)));
    hh:=h+a;
    a:=a/15;
    ra_dd:=ra_d-a;
    dec_dd:=arctan(cos(hh)*(r*sin(dec_d)-p_sin_l)/r*cos(dec_d)*cos(h)-p_cos_l);
    a_ra:=decimal_to_ra(ra_dd);
    a_dec:=decimal_to_degrees(dec_dd);
end;

    { ------------------------
        find coordinates of moon
        ------------------------ }

procedure find_crater(
            period:date_time;
            geo_lat,
            m_geo_long,
            m_geo_lat:dms;
    var m_long,
            m_lat,
            pos_angle:dms);
var
    m_geo_long_d,m_geo_lat_d,geo_lat_d,an_long,jd,t,bc,lc,c1,c2,o,f,a,
    x,y,sx,sy:double;
begin
    jd:=get_julian(period);
    t:=(jd-2451545)/36525;
    an_long:=125.044522-1934.136261*t;
    while an_long>360 do
        an_long:=an_long-360;
    while an_long<0 do
        an_long:=an_long+360;
    f:=93.271910+483202.0175*t;
    while f>360 do
        f:=f-360;
    while f<0 do
        f:=f+360;
    geo_lat_d:=degrees_to_decimal(geo_lat);
    m_geo_long_d:=degrees_to_decimal(m_geo_long);
    m_geo_lat_d:=degrees_to_decimal(m_geo_lat);
    bc:=arcsin(-cos(geo_lat_d)*sin(m_geo_lat_d)+sin(geo_lat_d)*
        cos(m_geo_lat_d)*sin(an_long-m_geo_long_d));
    y:=(-sin(m_geo_lat_d)*sin(geo_lat_d)-cos(m_geo_lat_d)*
        cos(geo_lat_d)*sin(an_long-m_geo_long_d));
    sy:=y;
    x:=cos(m_geo_lat_d)*cos(an_long-m_geo_long_d);
    sx:=x;
    a:=arctan(y/x);
    match_xy_with_degree(sx,sy,a);
    lc:=f-a;
    if lc>10 then
        lc:=lc-360;
    if lc<-10 then
        lc:=lc+360;
    c1:=arctan((cos(an_long-m_geo_long_d)-sin(geo_lat_d)/
        (cos(m_geo_lat_d*cos(geo_lat_d)+sin(m_geo_lat_d)*sin(geo_lat_d)*
         sin(an_long-m_geo_long_d)))));
    o:=ecliptic_obliquity;
    c2:=arctan((sin(o)*cos(m_geo_long_d))/(sin(o)*sin(m_geo_lat_d)*
        sin(m_geo_long_d)-cos(o)*cos(m_geo_lat_d)));
    m_long:=decimal_to_degrees(lc);
    m_lat:=decimal_to_degrees(bc);
    pos_angle:=decimal_to_degrees(c1+c2);
end;

end.